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ABSTRACT 

Radiation pressure on dust grains may be an important mechanism in driving winds 
in a wide variety of astrophysical systems. However, the efficiency of the coupling 
between the radiation field and the dusty gas is poorly understood in environments 
characterized by high optical depths like those in ultra-luminous infrared galaxies 
(ULIRGs) and massive dense star clusters. We present a series of idealized numerical 
experiments, performed with the radiation-hydrodynamic code ORION, in which we 
study the dynamics of such winds and quantify their properties. We find that, after 
wind acceleration begins, radiation Ray leigh- Taylor instability forces the gas into a 
configuration that reduces the rate of momentum transfer from the radiation field to 
the gas by a factor ~ 10 — 100 compared to an estimate based on the optical depth 
at the base of the atmosphere; instead, the rate of momentum transfer from a driving 
radiation field of luminosity L to the gas is roughly L/c multiplied by half the optical 
depth at the dust photosphere, which is far smaller than the optical depth in the deep 
interior. When we apply our results to conditions appropriate to ULIRGs and star 
clusters, we find that the asymptotic wind momentum flux from such objects should 
not significantly exceed that carried by the direct radiation field, L/c. This result 
constrains the expected mass loss rates from systems that exceed the Eddington limit 
to be of order the so-called "single-scattering" limit, and not significantly higher. 
We present an approximate fitting formula for the rate of momentum transfer from 
radiation to dusty gas through which it passes, which is suitable for implementation 
in sub-grid models of galaxy formation. Finally, we provide a first map of the column 
density distribution of gas in a radiatively-driven wind as a function of velocity, and 
velocity dispersion. 

Key words: galaxies: ISM — galaxies: star clusters — hydrodynamics — instabilities 
— ISM: jets and outflows — radiative transfer 



1 INTRODUCTION 



Dusty winds are ubiquitous in astrophysics: they are driven 
on scales ranging from single stars (e.g. 
star clusters 



(e.g. 



Habing||1996t to 

Lopez et al.||2011 ) to entire galaxies (e.g. 



Veilleux et al.| 20*05 ) . The driving mechanisms of these winds 
are diverse and in some cases uncertain, but one possible 
mechanism for many of them is the force exerted by radia- 
tion interacting with dusty matter. Photons moving through 
dusty gas can be scattered or absorbed by dust grains, trans- 
ferring some of their momentum. The grains, in turn, trans- 
fer this momentum to the gas either through hydrodynamic 
drag or via magnetic fields, possibly giving rise to a wind. 
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Radiation pressure on dust has been suggested as an impor- 
tant feedback mechanism in regulating star formation on 
the scales of galaxies ( Scoville||2003 Thompson et al.||2005 



[Andrews Thompson 2011) and individual massive stars 
clusters (O'del l et al.|1967||Scoville et al.|2 001; Krumholz 
Matzner 2009,jFall et al.|2010| [Murray et al.,2010, ,Krumholz 



& Dekel [20TO ), for driving dusty fountain flows in normal 



spirals (Chiao & Wickramasinghe 1972 



[ray et aL [[20051 [2011 [ [Hopkins et aL|2012 ) 



Elmegreen 


1983 


superwinds 1 


Mur- 



However, assessing these claims has been diflficult due 
to limited understanding of the radiation-matter interaction 
that drives the flow. For optically thin flows the problem is 
relatively simple, since the state of the radiation fleld is de- 
coupled from the gas. For optically thick media, however, 
the problem is signiflcantly more diflficult, because the gas 
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is capable of reshaping the radiation field. This can lead to 
a number of complex phenomena, such as photon bubbles 
(Blaes & Socrates 2003) and radiation Rayleigh- Taylor in- 
stability ( Krumholz et al. 2009 Jacquet K rumhol z||2011 



results and provide some caveats, and Section [6] summarizes 
our conclusions. 



Jiang et al.|2013 ). The existence of these behaviors makes it 2 MODEL SYSTEM 



non-trivial to calculate from first principles whether an ob- 
ject with a given set of properties can produce a radiatively- 
driven dusty wind at all, and, if it does, what properties that 
wind is likely to possess. This problem has thus far prevented 
definitive identification of the driving mechanisms for winds 
observed in a variety of systems (e.g. 
Quataert|2012| [Newman et al.||2012 ). 



Faucher-Giguere & 



In Krumholz & Thompson (2012 hereafter Paper I) we 
addressed the first part of this problem: under what circum- 
stances do we expect an object to launch a radiatively- driven 
dusty wind? We developed an idealized model system that 
allowed us to extract the important dimensionless numbers 
governing wind launching, and we then conducted numerical 
experiments with the radiation-hydrodynamics code ORION 
to explore the non-linear behavior of the system. The major 
results of Paper I are that the behavior of gravitationally- 
confined, dusty columns of matter subjected to radiative 
fluxes are governed primarily by two characteristic values: 
the dust optical depth and the Eddington ratio, both com- 
puted using the opacities that apply at the surface (i.e. the 
photosphere) of the dusty gas. For high optical depths and 
surface Eddington ratios close to but below unity, which may 
describe many galaxies and star clusters, we showed that 
radiation passing through the gas drives statistically steady 
turbulence with average Eddington ratio of unity, but not a 
wind. 

This result, however, does not answer the question of 
what happens if a wind is launched - either because the sur- 
face Eddington ratio exceeds unity, or because some other 
mechanism is able to eject matter, by itself or in conjunc- 
tion with radiation forces. This question is the main focus 
of our paper. We seek to determine at what rate the matter 
in a radiatively-driven wind is able to extract momentum 
from the radiation field, and how this depends on properties 
such as the strength of the radiative driving and the optical 
depth of the matter. In addition to illuminating the physics 
of the winds, we also derive a rough fitting formula that 
can be used in numerical simulations that do not include 
radiation- hydrodynamics, and instead treat radiative driv- 
ing using sub-grid semi-analytic models (e.g. Oppenheimer 
|fc Dave„2006l [Hopkins et al.^ 20TT| . piang et al [ ( |2013t per- 
formed preliminary work on this problem in the context 
of winds where the dominant opacity source is Thompson 
scattering from free electrons, and concluded that radiation 
Rayleigh- Taylor instability would limit the wind mass and 
momentum flux. We seek to investigate whether the same is 
true for dusty winds, and to extend their results by draw- 
ing quantitative rather than qualitative conclusions about 
how the wind momentum depends on the properties of the 
system. 

The remainder of this paper is as follows. In Section 
|2] we briefly review the basic equations and model system 
developed in Paper I, and consider how to extend them to 
the case of a dusty wind. In Sectionj3]we describe our nu- 
merical simulations, and in Section ^Jwe analyze the results 
they produce. In Section [5] we discuss the implications of our 



2.1 Governing Equations and Model System 

As in Paper I, we treat a section of a galactic disk or a 
young star cluster as an idealized model system consisting 
of a slab of gas with total surface density H filling the domain 
z > 0. A vertical radiation flux F — Fqz enters the domain 
of interest at z = 0, and there are no radiation sources at 
z > other than the thermal emission of the gas. The slab 
of material is confined by a constant vertical gravitational 
force per unit mass —gz] we neglect the self- gravity of the 
gas. 

Since we are interested in cases where the gas layer 
is optically thick, we describe this system using the two- 
temperature flux-limited diffusion (2TFLD) approximation, 
in which we assume that the radiation spectrum is locally a 
Planck function at every point, but we do not require that 
the temperature Tr describing this Planck function be iden- 
tical to the gas temperature Tg. In this approximation, in- 
teraction of radiation and matter is governed by the Planck 
and Rosseland mean opacities kp and kr. 



The equations governing this system are ( [Krumholz 
eraL][2007| 
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where v is the gas velocity, P — pksTg/ pmu is the gas pres- 
sure, p is the mean molecular weight in hydrogen masses, 
e = P/[{j — l)p] + i;^/2 is the gas specific energy, p is the 
mean mass per gas particle in units of the hydrogen mass 
mn, 7 is the gas ratio of specific heats, B = caTg /Atv is the 
frequency-integrated Planck function, E — aT^ is the radi- 
ation energy density, F is the radiation flux, A is the flux 
limiter, and R2 is the Eddington factor. In the 2TFLD ap- 
proximation, adopting the flux-limiter of [Levermore Pom-[ 
raning ( 1981|) and Levermore (1984), the radiation quanti- 
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= A + X'R. 

In this work we adopt opacities 

T 



(10 



-3/2 



10 



10 K 



cm g 



(8) 



(9) 



an approximation to the behavior of dust opacity at tem- 
peratures below ^ 150 K ( Semenov et al.|2003) . 

In Paper I we show that these equations are fuhy char- 
acterized by four non-dimensional parameters 



/e,* = 



Here kr > 



(10) 

(11) 

(12) 
(13) 



and similarly for A€p,*, and the 




starred quantities in turn are defined by 



(-) 

\ca J 



(14) 



We may think of T* as the characteristic temperature at the 
photosphere of the dusty gas, where the radiation escapes 
to infinity. The natural units of velocity, length, and time 
for the problem are 



9 



(15) 



In real astrophysical systems, 13 s is always very small for 
any non-relativistic flow, and ko is always of order unity 
and probably varies little from one galaxy to another. Thus 
in practice the quantities /e,* and r* determine the behavior 
of the system. 

2.2 Dimensionless Numbers for Winds 

If gas has been launched into a wind, it has obviously over- 
come its initial gravitational confinement. As discussed in 
the Introduction, in Paper I we show that this occurs only 
for /e,* > 1, i.e. only if the Eddington ratio at the dust pho- 
tosphere exceeds unity. In principle one wishes to determine 
the properties of winds launched at a range of /e,* > 1. 
However, we focus on the asymptotic limit /e,* oo, corre- 
sponding to a freely-accelerating wind with negligible gravi- 
tational confinement. Our reasons for doing so are threefold. 
First, this reduces the parameter space we must explore. 
Second, it is very likely that the case where there is no grav- 
itational confinement will produce the largest possible wind 
momentum flux, and so it can serve as a useful upper limit. 
We will see below that even this upper limit is quite restric- 
tive on the possible momentum of the wind. Third, in the 
case of radiating optically thick disks (e.g. starburst galaxy 
disks or the disks around QSOs), /e,* rises with height above 
the disk, so winds at large distances will have larger /e,* val- 
ues dZhang &: Thompson||2012| ). 

For /e,* oc, or equivalently ^ ^ 0, the quantities p*, 
/i*, and t* cease to be well-defined. It is therefore helpful 
to define alternative natural units in the freely-accelerating 
wind case. The sound speed Cs,* remains the natural unit 
of velocity, and to define a unit of time it is helpful to ask 



Table 1. Simulation Physical Paramters 



Name 


r* 


E 




ta 


/la/10-2 


Pa/10-16 






(g cm-2) 


(kyr) 


(kyr) 


(pc) 


(g cm-3) 


T3 


3 


1.4 


1.1 


6.9 


0.38 


1.2 


TIG 


10 


4.6 


1.1 


23 


1.3 


1.2 


T30 


30 


14 


1.1 


69 


3.8 


1.2 



Note that TIO describes both runs TIOLR and TIOHR, which 
have identical physical parameters but different resolutions and 
box size. All models have T* = 82 K, Cg,* = 0.54 km s--*^. 

Table 2. Simulation Numerical Parameters 



Name 


IC 




Lx X Lz 


Ax 


^run 


T3 


T3F0.5 


1024 X 16384 


85.3 X 1365 


0.083 


28.0 


TIOLR 


T10F0.5 


512 X 32768 


25.6 X 407 


0.05 


21.9 


TIOHR 


T10F0.5 


1024 X 16384 


25.6 X 407 


0.025 


10.6 


T30 




1024 X 16384 


27.3 X 437 


0.027 


3.6 



The initial condition (IC) column gives the name of the corre- 
sponding run in Paper I used to produce the initial condition. 
Nx X Nz is the size of the computational domain in cells, Lx x Lz 
is the size of the computational domain, Ax is the size of a com- 
putational cell, and Uun is the duration for which we run the 
simulation. All quantities are given in units of ha and ta- Note 
that all models were run during the phase where gravity was 
turned on with Ax/h^ = 0.5. Finally, for a description of the 
initial conditions for run T30, see the Appendix. 

how long it would take the momentum carried by the direct 
radiation field to accelerate matter from rest to this speed. 
The momentum flux per unit mass of the injected radiation 
field is 

Fo 



/rad,dir — 



Ec' 



and so we define the acceleration time as 

Cs , * Cs , * 



ta 



-U. 



(16) 



(17) 



/rad,dir K,R,*Fo/c /e,* 

If radiative trapping is significant, we expect the matter to 
increase its velocity by Cs,* in a time significantly shorter 
than ta . Finally, we can define characteristic length and den- 
sity scales from the combination of ta and Cs,*. These are 

h — + — ^* h — ^ — •^^'* 

/E,* ha r* 

We report all results in this paper in units of pa, ha, and ta 



(18) 



3 NUMERICAL SIMULATIONS 

We solve Equations ([T]) - Q with ^ = using the 
radiation-hydrodynamics code ORION. Our simulations are 
two-dimensional, and take place in the (x, z) plane; a flux F 
of radiation is injected at the bottom of the computational 
box, z = 0. The boundary conditions are periodic in the x 
direction, impermeable at the lower z boundary and open at 
the upper z boundary. More details on the boundary con- 
ditions are given in Paper I. All other parameters of the 
simulations are also the same as in Paper I, except that we 
set the external gravitational field ^ = 0, so that /e,* oo. 
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To study the behavior of winds in the hmit /e,* oc 
as a function of r* , we run four simulations, which we denote 
T3, TIOLR, TIOHR, and T30; the first of these have r* = 3, 
the second two have r* = 10, and the third has r* =30. The 
two T* = 10 runs are identical except in their resolutions and 
the sizes of the computational domains. We perform both 
runs as a check on the resolution-dependence of our results. 
We summarize the physical and numerical parameters of the 
simulations in Tables [T] and [21 

In Paper I, we show that initially laminar, 
gravit at ionally- confined gas slabs subjected to radia- 
tion fluxes such that /e,* is above a certain critical 
value, but less than unity, will develop radiation-driven 
turbulence. This occurs in a time < lOOt*. (Note that 
values of t* for our runs are given in Table [l] and are the 
same for every run; lOOt* ^ 100 kyr.) Since this is short 
compared to most astrophysical time scales of relevance, 
it is reasonable to assume that gas being launched in a 
wind will be in a fully turbulent state. We therefore do not 
start our simulations with laminar gas layers. Instead, we 
use as initial conditions the end states of the simulations 
from Paper I, as summarized in Table [2] We modify these 
conditions only in that we place the gas in a computational 
box that is larger in the vertical direction, in order to 
accommodate vertical expansion of the gas layer once it is 
no longer gravitationally-confined. We initialize computa- 
tional cells that are outside the computational domain of 
the simulations of Paper I by giving them densities equal 
to the background density and temperature from Paper 
I, 10~^°y9* and T*, respectively. In run TIOLR, we also 
down-sample the resolution by a factor of 2. The exception 
to the above statements is run T30, for which we do not 
have a corresponding run from Paper I. We describe how 
we generate its initial conditions in the Appendix. 

As the simulations proceed, when necessary we shift 
all velocities in the computational domain by a constant 
offset in order to bring the center of mass velocity of the 
gas back to zero. Our method is simple: we have added an 
option to the ORION code that, upon restart from a check- 
point, calculates the center of mass velocity of the compu- 
tational domain in the z direction, then subtracts the corre- 
sponding velocity from all computational cells, altering the 
momenta and total energies appropriately. The calculation 
then restarts from the modified data. We apply this option 
whenever a visual inspection of the data indicates that the 
bulk of the mass is well away from the bottom boundary of 
the computational domain. This enables us to continue the 
simulations longer without the gas reaching the top of the 
computational box. We shift the velocities in this manner 
only when the vast majority of the gas is well away from the 
bottom of the computational box, so that there are no sig- 
nificant artificial forces exerted by the bottom of the compu- 
tational box. In the analysis below, we remove these offsets 
and present the results as if the entire simulation had simply 
taken place in a larger box. In principle we could shift the 
positions as well, but this is less convenient computation- 
ally, since it would require translating values from one cell 
to another, and filling in values of density, momentum, and 
energy in the new cells added to the computational domain 
by any shifts. 



4 RESULTS 

4.1 Qualitative Behavior 

Figures [l]-|8] shows a series of snapshots of the simulation 
density fields. As the plots show, the initial state that re- 
sults from the radiation Ray leigh- Taylor (RRT) instability 
( [Jacquet k, Krumholz 20TT| acting on a gas confined by 
gravity, consists of a relatively horizontal, turbulent layer. 
In the absence of gravitational confinement, the radiation 
force rapidly drives the gas into a predominantly vertical, 
filamentary structure. In between the filaments of dense gas 
there are low-density channels. As the material is accelerated 
upward by the radiation field, the gas becomes more elon- 
gated and spread over a progressively larger vertical extent. 
We are eventually forced to halt our simulations primarily 
because the vertical extent of the gas becomes comparable 
to the vertical size of our computational domain. 

Figure [9] shows an example of the distribution of density, 
temperature, velocity, and radiation flux in one of the runs 
once the channel structure has developed. The channels are 
characterized by several features. First, within them the gas 
is traveling at extremely large velocities relative the dense 
gas in the filaments. At the snapshot shown, the velocity 
difference approaches many tens of Cg,*. As a result of this 
velocity difference, the edges of the channels appear to be 
scalloped by Kelvin-Helmholtz instabilities. Second, because 
of their lower optical depths, the channels carry the great 
majority of the radiative flux. The flux within the chan- 
nels approaches lOFo, while inside the filaments the flux is 
^ Fq. Thus the matter effectively collimates the radiation 
field, inducing a strong anti-correlation between density and 
radiative flux. This anti-correlation is the main signature of 
the RRT instability. 



4.2 Radiative Trapping 

The development of vertical filamentary structure and the 
resulting coUimation of the radiation field has profound ef- 
fects on its ability to trap the radiation field and extract 
momentum from it. To quantify the rate at which the gas 
takes up momentum from the radiation field, it is helpful 
to examine the z component of the momentum equation in- 
cluding radiation and gravitational forces; this is 



dP , Fz 
dz c 



P9, 



(19) 



where Fz is the z-component of the radiation flux and the 
use of kr in the equation implicitly equates the flux-mean 
and Rosseland-mean opacities, as is appropriate in the opti- 
cally thick regime. In the flux-limited diffusion approxima- 
tion this equation is equivalent to Equation ^ , as shown by 
iKrumholz et al. (2007), but the analysis is more transparent 



when the equation is written in the form above. If we inte- 
grate this equation over the entire computational domain, 
and ignore the small terms that arise from forces and fluxes 
across the top and bottom boundaries of the computational 
domain, the first two terms on the right-hand side vanish 
and we are left with 



dt 



{pVz 



I^RP- 



{p)9, 



(20) 



where for any quantity q we defined the volume average by 
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Figure 1. Density distribution as a function of time in run T3. 
Snapshots are shown at intervals of 5ta, starting from t = 0, as 
indicated at the top of each panel. White bars indicate a region 
around the vertical center of mass; we show a zoom-in of this 
region in Figure [s] Note that the vertical extent shown does not 
necessarily match the size of the computational box given in Ta- 
ble [2] because we have compensated for the effects of our periodic 
shifts the center of mass velocity of the entire computational do- 
main in some runs - see Section |3] for details. 



Figure 2. Same as Figure^ but for run TIOLR. The zoomed- 
region is shown in Figure [6] 
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Figure 3. Same as Figure^ but for run TIOHR. The zoomed-in Figure 4. Same as Figure^ but for run T30. The zoomed-in 

region is shown in Figure [7| region is shown in Figurejs] Note that the first panel is not in fact 

empty - the gas at time is simply compressed into an extremely 
thin layer whose width, on the scale plotted, is less than a single 
pixel. 
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-42 -42 -42 -42 -42 42 

x/K x/h^ x/h^ x/h^ x/h^ 



700 



600 



500 



Figure 5. Same as Figure [T] but the panels show a zoom-in on 
a smah region around the vertical center of mass. The zoomed 
region is indicated by the white bars in Figure^ 



400 



LxLz 



Jo 



q dz dx. 



(21) 



Dividing both sides by (p) and noting that the mass in the 
computational domain is very close to constant (since we 
are careful to ensure there is no significant mass loss from 
the top of the computational box) gives 



dvz 
dt 



1 (i^RpFz) 



-9, 



(22) 



(23) 



where Vz = {pVz)/{p) is the mass- weighted mean z velocity 
of the gas. We use the first term on the right-hand side to 
define the mean radiation force per unit mass, 

1 {krpFz 
' c ip) 

Equivalently, we may think of this term as describing the 
mass- weighted mean radiation force. Based on our observa- 
tion that density and flux are strongly anti-correlated, we 
expect that to be much less than the volume- weighted mean 
radiation force (krFo)/c would be. The second term on the 
right-hand side is simply the gravitational force per unit 
mass. 

At this point it is useful to rewrite the equation by mul- 
tiplying through by a factor of ta/cs,* to non-dimensionalize. 
Doing so gives 

ta_^^ _ /rad 

Cg * dt 



/rad,di] 



/e,* 



(24) 



Following Krumholz & Matzner (2009) and Krumholz & 
Thompson (2012), we define the trapping factor by 

l+/t.ap = y^. (25) 

/rad,dir 

Physically, the trapping factor is simply the factor by which 
the radiation force is amplified by trapping of the radiation 
field by the gas. The quantity 1 + /trap is equivalent to the 



300 



200 



100 



-12 



12 



Figure 6. Same as Figure [s] but for run TIOLR. 



amplification factor tir defined by Thompson et al. (2005), 



although we refer to it as /trap here because, as we will see, 
its relationship to optical depth is not trivial. Using equation 
(24) to rewrite the equation above, we obtain 



/trap — 



tg dVz 

Cs * dt 



/e,* 



(26) 



The quantity on the right-hand side is directly measurable 
from our simulations (and /e,* = oo in the absence of grav- 
ity, so the term t*//e,* = 0), so our simulations provide 
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t=0 t=2.5 t=5 t=7.5 t = 10 




-12 -12 -12 -12 -12 12 

x/K ^IK ^IK x/ha ^IK 



Figure 7. Same as Figure [s] but for run TIOHR. 



t=Q t = l.l7 t=2.33 t=3.5 




x/K x/K x/h^ x/h^ 



Figure 8. Same as Figure [s] but for run T30. As with Figure^ 
note that the first panel is not empty, but it appears so because 
the gas at time is simply compressed into an extremely thin 
layer whose width, on the scale plotted, is less than a single pixel. 




Figure 9. A section from run TIOLR t = 12Ata showing 
the density and velocity distribution (colors and vectors in the 
left panel), and the gas temperature and radiative fiux distribu- 
tion (colors and vectors in the right panel). The region shown is 
centered in the vertical center of mass of the gas at this time, 
Zcm = 181/ia, and the velocities shown are relative to the vertical 
center of mass velocity of the section shown, 

'^2, cm — 

77.0c...*. 



us with a direct measurement of /trap as a function of time. 
However, we must make one important modification to equa- 



tion (26), which comes from a limitation of our numerical 



method. Because we are using flux-limited diffusion, we do 
not properly capture the interaction of the gas with the di- 
rect, beamed radiation field produced by stars. Instead, we 
are treating the radiation field only after this first absorp- 
tion. Since the final —1 represents the contribution from 
this direct radiation field, we do not subtract it off when 
computing /trap from the simulations. This is likely conser- 
vative, since our method does capture some of the effects 
of the first absorption, in which case the results we obtain 
should be upper limits on /trap- However, we cannot com- 
pletely rule out the possibility that inclusion of the direct 
radiation force would somehow change the structure of the 
gas and indirectly increase the trapping of the reradiated 
field. 

Figure ^] shows the gas mean velocity and trapping 
factor as a function of time in each of our simulations. For 
constant /trap, the gas velocity should increase linearly with 
time. Instead, we see that the velocity increase is steep at 
first and then becomes much shallower, and this is reflected 
in the plots of /trap, which are large at first and then decline 
over a few ta- The initially high values are easy to under- 
stand given our starting conditions. When the gas is confined 



by gravity and there is no wind, both the time- averaged 
value of Vz and its rate of change must be zero. Consult- 
ing equation (26), this requires that /trap = t*//e,*, again 
omitting the —1 because our simulation does not properly 
model the direct radiation force. The initial value of /trap 
we measure is indeed close to this, though the match is not 
exact because dvz/dt is not precisely zero at all times in the 
gravity- confined state; instead, it oscillates about zero. 

Once the gravitational confinement is removed, how- 
ever, the gas morphology changes from predominantly hor- 
izontal to predominantly vertical, and /trap drops. At late 
times /trap oscillates up and down about a value well below 
the initial one. Upward and downward oscillations of /trap 
correspond to variations in the gas morphology. At times, for 
example at times t/ta — 5 and 10 in run TIOHR (see Figure 
|7|, the filaments formed by the radiation are fairly coher- 
ent leave fairly large vertical channels unobstructed, and at 
these times /trap is low. At other times, such as t/ta — 7.5 
in run TIOHR, the filaments are more fragmented and cover 
more or the domain horizontally, giving rise to larger values 

of /trap ■ 

Finally, comparing runs TIOLR and TIOHR suggests 
that are results are relatively well- converged. Since the gas 
is turbulent for RRT instability in its fully developed state, 
the flow is chaotic and we do not expect either morphologies 
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Table 3. Simulation Results 


Name 


^min 


min(/trap) 


(/trap) 


max(/trap) 


/ /e,* ^ /trap,grav 


l^{Tmp)^ ^ TiR 


T3 


5 


1.0 


1.6 ±0.4 


2.5 


6 


15 


TIOLR 


3 


2.6 


6.0 ± 1.3 


9.9 


20 


120 


TIOHR 


3 


2.2 


4.9 ±2.4 


8.9 


20 


120 


T30 


1 


5.3 


12.3 ±5.7 


23.3 


600 


2000 



For each run, min(/trap), (/trap), and max(/trap) give the minimum, time-averaged, and max- 
imum values of /trap that occur in the simulation after time tmin- The error bars given on 
(/trap) represent the la range measured from the simulations at times > tmin- For comparison, 
't^/Ie,* ~ /trap,grav givcs the timc-avcragcd value of /trap in the steady-state gravitationally- 
confined configuration from which we start, while ^t(Tmp)S ~ tjr is the average optical depth at 
the start of the calculation, computed using the mass-weighted mean midplane temperature as in 
Paper I. 




25 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 



Figure 10. Mean gas velocity versus time (top row) and trapping 
factor versus time (bottom row) for runs T3 (first column), TIO 
(second column), and T30 (third column). In all panels solid lines 
show the results of the simulations; for TIO, blue shows TIOLR, 
and green TIOHR. In the upper panels black dashed lines show 
the change in velocity versus time that would be expected for 
trapping factors /trap = 1, 10, 100, and 1000, from shallowest line 
to steepest. In the bottom panel, dashed black lines show /trap = 
1 and /trap = t*//e,*, with the value of /e,* computed before 
gravity is turned off - i.e. the value of /e,* that was used in the 
simulation from Paper I from which we take our initial conditions. 
In the absence of gravity, as is the case for the simulations shown 
here, /e,* = oo. 



0.10 
; 0.08 
I 0.06 

i 0.04 

I 

: 0.02 

0.00 

; 0.08 

i 0.06 
i 0.04 

I 

: 0.02 

0.00 

; 0.08 

I 0.06 
I 0.04 

: 0.02 

0.00 

; 0.08 

i 0.06 
i 0.04 

I 

: 0.02 



0.00 



-20 



-50 



20 



[km s M 
40 60 



80 100 120 



T3 




TIOLR 



■ TIOHR 



T30 



t/ta =0 
t/ta =3 



50 



100 



150 



200 



250 



or exact values of /trap as a function of time to be resolution- 
independent. However, examining the results in Figure 
we see that values of /trap versus time produced in the two 
runs are qualitatively similar, and that quantitatively their 
means are well within the level of variance in /trap we mea- 
sure in each run as the flow varies chaotically. This suggests 
that our values of /trap are converged. 

We summarize our results for /trap in Table [3] where 
we report the minimum, mean, and maximum values of 
/trap we measure in each of our simulations once the ini- 
tial transient phase ends. For comparison, we also report 
t*//e,* ~ /trap,grav, the mcau value of /trap in the initial, 
turbulent, gravity- confined state, and niTmp)^ ~ tir, the 
optical depth computed by multiplying the column density 
by the opacity evaluated using the midplane temperature 
Tmp. The latter has been used as an approximate value for 
/trap by a number of authors, as we discuss in more de- 



Figure 11. Velocity distribution functions for each of our sim- 
ulations at times t/ta = 0, 10, and 28 (for runs T3, TIOLR, and 
TIOHR) and t/ta =0 and 3 (for run T30), as indicated in the leg- 
end. In each panel the histogram shows the fraction of the mass 
in the simulation that falls into a given bin of z velocity at the 
indicated time. Vertical dashed lines indicate the mass-weighted 
mean velocity Vz at that time. 



tail below. Clearly none of these values are equal; instead 

/trap /trap,grav < TIR. 

4.3 Wind Velocity Distribution 

Figure [To] and Table |3] describe the mean velocity and mass- 
averaged momentum transfer from radiation to gas. How- 
ever, it is also interesting to look at the distribution of mat- 
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ter velocities. In Figure [TT] we show mass- weighted velocity 
probability distribution functions (PDFs) for each of our 
runs at several times. In the initial condition, the distribu- 
tion of z velocities is, as one might expect for a roughly hy- 
drostatic atmosphere, symmetric about Vz = The width 
is fairly small, reflecting the relatively low Mach numbers 
we obtained for turbulent atmospheres in Paper I. At times 
^ ta, not only is the mean velocity Vz larger, the spread 
of velocities is larger as well. At late times the PDF is all 
the runs is slightly asymmetric, with the majority of the 
mass residing at velocities slightly below the mean, and a 
tail extending well above the mean. 

The division between high and low velocity material 
corresponds to the division between material in the opaque 
filaments and material in or at the edges of the radiation- 
dominated channels, as illustrated in Figure[T2] To construct 
this figure, along every vertical line of sight we measure the 
column density 



E(x) 



/' 



z) dz. 



(27) 



We then assign every cell a column density correspond- 
ing to the value at its x position, and construct the two- 
dimensional PDF of and Vz. From the 2D PDF, we 
see the same asymmetry as in Figure [111 where the veloc- 
ity distribution extends further from the mean in the posi- 
tive direction than the negative direction. In the 2D PDF, 
it is clear that the high velocity material consists prefer- 
entially of gas with low The correlation is relatively 
weak, and the overall range in is relatively small, be- 
cause the filaments are not perfectly vertical. Thus, most 
of the time a given vertical line of sight will intersect both 
dense filaments and low-density channels, rather than look- 
ing straight down the barrel of a channel. Nonetheless, this 
column density-velocity anti-correlation represents a possi- 
ble observable signature of radiation pressure- driven dusty 
winds. 

The overall width of the velocity distribution, including 
both low and high speed components, is of order ^ 20cs,* 
in all the runs; the dispersion of horizontal velocities is sub- 
stantially smaller. The dispersion does not appear to in- 
crease substantially over the time interval shown, and thus 
^ 20cs,* is likely the steady-state value, at least over the 
range of r* values that we have explored. This corresponds 
to a one-dimensional Mach number in the vertical direc- 
tion of order 20 - not exactly 20, since much of the gas 
is somewhat warmer than T* and thus has a sound speed 
greater than Cs,*. For our fiducial choice of dimensional 
scaling (cs,* = 0.54 km s~^), this give a physical velocity 
dispersion of roughly 10 km s~^ in the wind, compared to 
bulk velocities of ^ 100 km s~^ at the same time. While 
this Mach number and velocity dispersion are larger than 
we found in Paper I for RRT-unstable atmospheres that do 
not drive a wind, they are still close to an order of magni- 
tude smaller than the values observed in the most vigorously 



star- forming ultraluminous infrared galaxies (e.g. |Downes 



-3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 

v^-v^ [km s"^ ] 
-30 -20 -10 10 20 30 40 




Solomon 1998). Finally, we note that the results for runs 
TIOLR and TIOHR are qualitatively similar at equal times, 
suggesting at least rough convergence. 



-20 20 



Figure 12. Two dimensional velocity-column density distribu- 
tions in the simulations at the latest time slices shown in Figure 
[Tl](V^a = 28 for runs T3 and TIOLR, t/ta = 10 for run TIOHR, 
and t/ta = 3 for run T30). Each pixel shows the logarithm of 
the fraction of the simulation mass in the indicated bin of Vz 
and ^(x), normalized so that the most massive bin has a value 
of unity. Note that both the x and y axes are offset such that 
material at the mass-weighted mean column density and velocity 
would appear at the position (0,0). 



5 DISCUSSION 

5.1 Fitting Formulae for Radiation Trapping 

By combining the results of this paper with those of Paper I, 
we are now in a position to provide a fitting formula for the 
value of /trap in optically-thick radiation pressure-driven at- 
mospheres and winds. Such a formula is useful in simulations 
or analytic calculations that seek to include radiation pres- 
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Figure 13. Values of /trap as a function of r* measured in 
the simulations presented in this paper. The points represent the 
time-averaged value, thick error bars show the standard deviation, 
and thin error bars show the range from minimum to maximum; 
all values are as given in Table [3] For runs TIOLR and TIOHR, 
the points are offset slightly from r* = 10 for clarity. The dashed 
black line is /trap = 0.5r*. 



sure effects, but that do not properly capture the radiation- 
hydrodynamic behavior associated with the RRT instability. 
Examples include one-dimensional models (e.g.'Krumholz & 



[Matzner 2009; Murray et al. 2010), models that include ra- 
diation pressure only via a subgrid prescription rather than 
by solving the equation of radiative transfer (e.g. Hopkins 
et al.|20lT Agertz et al.|2012 ), models that solve the trans- 



fer equation only for the direct and not the dust-reprocessed 
radiation field (and therefore implicitly set /trap = 0; e.g. Pe- 
ters et aLjlWo] [Wise et al.||2012| |Kim et al.||2013a|bt , and 



models that solve the transfer equation in one dimension 
under an assumption of spherical symmetry and therefore 
miss RRT effects (e.g. [Novak et al.|2012 ). 

The value of /trap for a radiation pressure-dominated 
wind or atmosphere is a function of the two main dimen- 
sionless parameters for the problem, /e,* and r*. We have 
sampled this parameter space quite coarsely, but we can 
nonetheless provide a rough fit that captures the results of 
our simulations, and which is an improvement over simple 
prescriptions. In Paper I we explored the regime /e,* < 1, 
and found that for a given r* there exists a critical /e,* above 
which instability sets in. At values of /e,* below this value 
the atmosphere is supported predominantly by gas pressure 
over most of its height, and radiation pressure is dynami- 
cally unimportant. Above the critical /e,* values, we found 
that RRT instability causes the value of /trap to self- adjust 
so that the radiation force exactly balances gravity without 
producing a wind. This is /trap ~ t*//e,* — 1. To extend 
this to lower r* than we have sampled, we simply impose 
the requirement that /trap cannot be less than 0. Thus our 
approximation for /e,* < 1 is 



/trap,io ~ max 



1,0 



(28) 



In the regime /e,* oo that we explore in this paper 
there is a wind, but the rate at which it takes up momen- 




-0.5 0.0 0.5 

log /e,* 



Figure 14. Values of log /trap (top), log{/E) (middle), and 
{dPwind / dt) / {L / c) (bottom) as a function of /e,* and r*, com- 
puted using the fitting formula given by Equation ( |30| . In each 
panel contours lines appear at values of —2, —1, 0, 1, and 2. The 
thick black line in the top two panels shows the critical curve 
below which RRT instability shuts off (white region). The thick 
black line in the bottom panel shows the critical value below 
which no wind is launched. 



tum from the radiation field is limited. Figure [13] shows our 
estimated values of /trap as a function of r*, together with 
a crude linear fit that is consistent with the simulations: 



/trap, hi ~ 0.5t* 



(29) 



This is only a "by-eye" fit, but it describes the data very well, 
and given the small number of simulations and the error bars 
on each one, a more sophisticated fitting procedure does not 
seem justified. 

To combine the two cases /e,* < 1 and /e,* oo, 



we hypothesize that /trap will obey equation ( 28 ) up to the 
point where /e,* = 1. At this point we will have /trap ~ r* — 1 
for T* ^ 1. Beyond this point, as /e,* increases /trap will 



smoothly decrease onto the fit given by equation ( 29 ) in the 
limit /e,* oo. Since we have not mapped out the inter- 
mediate /e,* regime, obviously the functional form of this 
decrease is not well- constrained, and we cannot rule out the 
possibility that /trap behaves non-monotonically over this 
range, for example developing a peak at a special value of 
/e,*. However, there is no good reason to believe that such 
a phenomenon should occur, and in its absence the func- 
tional form we adopt to interpolate between the behavior 
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at /e,* < 1 and /e,* oo matters little, since /trap only 
changes by a factor of two between those cases. We therefore 
adopt a very simple interpolation between the two cases, 

/trap,lo /trap, hi 



/trap ^ /trap, hi H~ 



max(/E,*, 1) 



(30) 



in the regime where RRT instability occurs. This fitting 
formula has the correct asymptotic behavior in the limits 
T* ^ 0, /e,* < 1, and /e,* oo, and is consistent with the 
simulations over the range of /e,* and r* we have explored. 
We can also use this formula to determine the ratio of mass- 
averaged radiation force to gravitational force. In Paper I 
we showed that this is 



(/b) = (1 + /trap) 



(31) 



If a wind is launched, from equation (24), we can see that 



the net rate at which it gains momentum including both 
radiative acceleration and gravitational deceleration is 



wind _/-,,/> \ / (/e 
- (1 + /trap) 



(/b> 



(32) 



Figure 



14 



shows the results of our fit for /trap, (/e), and 
(dp^ind / dt) / {L / c) as a function of r* and /e,*- In the plot 
we can see the three dominant regimes identified in our sim- 
ulations. Below the critical line there is no RRT intst ability, 
and radiation is dynamically subdominant. For /e,* large 
enough to turn on RRTI but still below about unity, (/e) is 
fixed to unity, and /trap self-adjusts to compensate, decreas- 
ing as /e,* increases at fixed r*. Finally, at /e,* greater than 
about unity, a wind appears. In this regime /trap is a func- 
tion primarily of r*, and is quite insensitive to /e,*- On the 
other hand (/e) increases with /e,*, indicating that gravity 
is becoming progressively weaker relative to radiation. As a 
result, the wind strength is monotonically increasing with 
/e,*, but only slowly, since gravity is relatively unimportant 
once /e,* is even a slightly above unity. We stress that the 
exact location of the wind-launching line should not be taken 
too seriously, particularly at r* < 1, given the sparsity with 
which we have sampled the parameter space. Nonetheless, 
the qualitative result that for r* > 1 a wind appears only 
for /e,* ^ 1 should be robust. 

Readers may note that, for r* <^ 1, it is possible for 
there to be a wind even when /e,* < 1 and RRTI does 
not occur. Physically, this corresponds to a medium that is 
optically thin to dust-reprocessed radiation, but is still ab- 
sorbs the direct radiation field. In this case, the radiation 
is absorbed once, is reemitted, and then immediately es- 
capes, so the RRTI that we see in our simulations does not 
occur. However, if the direct radiation field carries enough 
momentum, this single absorption may still be sufficient to 
overcome gravity and launch a wind. Simulations by (Kuiper| 



et al. (2012) suggest that RRTI does not occur in this case, 



which is not surprising, since RRTI relies upon the ability 
of the gas to shape the radiation field. That cannot happen 
if the radiation is only absorbed once. 



5.2 Limitations Due to Geometric Simplifications 
of the Simulations 

Our simulations represent an idealized numerical experiment 
with a simple geometry. It is therefore interesting to ask how 



a more realistic setup would likely affect our results. One 
obvious simplification in our simulations is that they are 
two- rather than three-dimensional. The implications of this 
are discussed extensively in Paper I, and we refer readers to 
the discussion there. 

A second simplification is that we have assumed a pla- 
nar geometry, whereas a real wind will generally approach 
a spherical geometry, at least once it is far from its launch 
point. A small section of a spherical shell of wind material 
may be treated as planar, and so our planar results should 
continue apply locally. The main difference between planar 
and spherical geometries, therefore, is that in planar geom- 
etry T* and /e,* are fixed, whereas for a spherical wind they 
will vary as the wind expands. This variation is caused by 
two effects. First, as a spherical shell of constant mass ex- 
pands in radius i?, its surface density drops as R~^, which 
reduces r*. Second, the gravitational force g and the flux Fq 
encountering a spherical shell also both drop as . This 
means that the ratio Fo/g remains constant; however, the 
drop in Fq reduces T* as R~^^'^ and thus kr^^ as which 
affects both /e,* and r* . The net effect is that, for a spherical 
shell of fixed mass and radius i?, /e,* oc and r* oc R~^ , 
and thus an expanding spherical shell traces a line of slope 3 
in the (log/E,= 



,logr*) plane depicted in Figure 14 Systems 
start at the upper right of the plane, then move down and 
to the left as they expand. Examining the Figure, we see 
that such a trajectory will result in a value of /trap that de- 
creases with time, and that it approaches an asymptotically 
constant value of (/e)- This is not surprising; it is simply a 
statement that, as a shell expands, its optical depth drops 
and thus the dust-reprocessed radiation field becomes less 
and less important compared to the direct one, which gives 
constant (/e)- 

A third simplification of our simulations is that we have 
assumed a constant flux as would would be produced in 
a galaxy with all the stars at the midplane, or by a star 
cluster with all the stars concentrated at the center and the 
mass at fixed radius. In reality, the sources of radiation are 
intermixed with the gas being launched in a wind. This may 
well result in a significant reduction of the direct radiation 
force due to geometric cancellations. However, it should not 
substantially affect our results for /trap, simply because an 
appreciable value of /trap requires that the radiation field 
be trapped and therefore isotropized, forgetting its original 
direction. Thus our results for /trap should be robust against 
a change from planar or point-like sources to distributed 
sources. 



5.3 Implications for Star-Forming Systems 

Given our results for /trap and (/e), it is interesting to ask 
what our models predict for star-forming systems, which 
have been posited to be regulated by radiation pressure. We 
consider two types of objects: proto-star clusters, which we 
approximate as spherical, and galactic disks, which we ap- 
proximate as planar. For a spherical object of total (gas plus 
stellar) mass M, gas mass fraction /^, stellar mass fraction 
/* = ^ — fg, and radius R, within which the stars have a light 
to mass ratio ^, we have a central luminosity L — "^f^M, 
gas surface density S = (1 — f^)M / A-kR? , surface gravita- 
tional acceleration g — GM/R^ , surface flux F — L/AttR^, 
and surface temperature T* = (L/47ri?^cr)"^/^. Thus for the 
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rough fit to the [Semenov et aL ( |2003 ) opacity model given 
by equation ([9|, our estimates of the key dimensionless pa- 
rameters near the surface become 



r* = 2.6/^,0.5/ 



1/2 



0.5^3 



1/2^3/2 



E,* 



3/2^1/2 



(33) 
(34) 



where /*,o.5 = /*/0.5, /,,o.5 = (1 - /*)/0.5, Si = S/10 g 
cm~^, and ^3 = ^/IO^Lq/M©; for a zero-age stellar pop- 
ulation, ^ = II4OL0/M0 f Fall et al.||2010 '), so our nor- 
malization should be appropriate for a young cluster. In 
computing r*, we have assumed that the gas is arranged 
in a thin shell (as it must be if it is to be ejected), so 
T* = i^^M / A^nR^ — Ai:*S/4, with the factor of 4 arising from 
the difference between the surface density S as measured by 
an external observer (which is what appears in the above 
equations) and the surface density as seen by radiation es- 
caping from the center of the star cluster; if we instead adopt 
a uniform density sphere geometry, r* will be larger by a fac- 
tor of 3. Similarly, variations in the dust opacity per unit gas 
mass could plausibly increase r* and /e,* by as much as a 
factor of a few, and downward by much larger factors in low- 
metallicity systems. Note that, as pointed out by [Fall et al.| 
( 2010 ), M and R enter only through the combination E. The 
value to which we have normalized S, 10 g cm~^, is roughly 
the maximum observed value for stellar systems anywhere 
in the Universe ( Hopkins et al.|20'l0 ), and thus the values of 
r* and /e,* above should be regarded as upper limits. Note 
further that, even though /e,* < 1 in equation (34), this 



does not imply that radiation is unimportant to the dynam- 
ics of the dusty gas. Indeed, radiation pressure should drive 
strong turbulence as in our simulations presented in Paper 
I. 

We can perform a similar calculation for a galactic disk 
with total surface density S and gas mass fraction fg. Such 
as disk has a surface gravitational force g — 27iG^. For 
a stellar population older than ^ 4 Myr, the light to star 
formation rate ratio approaches a roughly constant value 
$ = 6.1 X 10^^ erg g"^ s"^ = 1.0 x 10^° L0/(Mq yr"^) (cal- 
culated using starburst99 -|Leitherer et al.|1999| |KrumhoIz| 
k Tan||2007| [Krumholz Dekel||2010t . It is therefore con- 



venient to write the radiative flux as F = where 
is the star formation rate per unit area. Plugging this flux 
into our scalings for the dimensionless parameters gives 



= 0.67/s,o.5$l^' 

3/2^3/2 



o.43$J/^s: 



3 ^0 



(35) 
(36) 



where So = ^/l g cm"^ S*,3 = /lO^ Mq pc"^ Myr"\ 
and $10 = $/10^° Lq/Mq. The normalizations of S and 
here have been chosen to match those of the most vig- 
orously star-forming galaxies observed. Indeed, none of the 
galaxies in the large sample compiled by [Krumholz et al!] 
(2012) exceed this star formation rate. Since observations 
indicate that oc with p 1 — 1.5, using a normaliza- 
tion for galaxies of lower star formation rates and surface 
densities would lead to lower values of r* and /e,*- Thus 
the values above are, as in the case of single clusters, upper 
limits for star-forming systems. However, we note that QSO 
disks on ~ 1 — 50 pc scales can and do exceed these limits 
( Sirko Goodman||2003| [Thompson et al.||2005 ). 

Given these numbers for star clusters and galactic disks, 
we can draw a few conclusions. The first, already suggested 



in Paper I, is that, in the absence of additional forces, the 
dust-reprocessed radiation field cannot launch winds or eject 
mass from the great majority of star clusters and galaxies. 
This is because we find that winds are only launched when 
/e,* > 1. For star clusters even our upper limit is well below 
this value, and for galaxies only the most extreme systems 
approach it, while galactic winds are inferred to be ubiq- 
uitous (e.g. Veilleux et al. 2005). This is not to say that 



radiation pressure is not important. As discussed above, if 
the radiation force is sufficiently strong and r* < 1, it may 
be possible for the direct radiation field to eject matter, par- 
ticularly in systems where gravity is already partially offset 
by magnetic fields, turbulent motions, or some other force 



(e.g. [Murray et al.||2005| jKrumholz k Matzner||2009| [Mur- 



ray et al.|2010||Fairet al.|2010| [Hopkins et al.|2011t . Indeed, 
Krumholz & Matzner ( 2009 ) compile a sample of super-star 



clusters, and show that for some of them the direct radiation 
force, combined with the momentum of line-driven winds, is 
likely to be able to eject matter even without significant 
radiative trapping. Even in somewhat lower luminosity sys- 
tems where it cannot eject the bulk of the gas, radiation 
pressure may still be able to drive small amounts of mass to 
speeds above the escape speed and eject it, as happens for 
example with massive stars. Nonetheless, our results show 
that ejecting mass from star-forming systems via radiation 
pressure is significantly more difficult than many models as- 
sume. 

A second implication of our work is that, if radiation 
pressure does launch winds, and if it were the sole driv- 
ing mechanism, those winds are not likely to carry a mo- 
mentum flux much larger than a few times L/c. Equations 
(33) and (35) show that r* on galactic scales never much 



exceeds unity, and that even for the densest clusters it is 
< 10; a more typical value for massive clusters would give 
r* ^ 1. Since we find that /trap ~ 0.5t*, this means that we 
cannot expect winds accelerated primarily by radiation to 
have /trap larger than ^ 1. Thus radiation pressure-driven 
winds from star-forming systems should not carry a mo- 
mentum flux that exceeds L/c by more than a few tens of 
percent. The best fit values of the momentum fluxes of the 
winds produced by giant star-forming clumps at z ~ 2 ex- 
ceed this limit (Genzel et al. 2011| Newman et al. 2012), 
which taken at face value would suggest that they cannot 



be primarily radiation-driven (consistent with Krumho lz fc] 
^Dekel 2010) . However, we caution that there are very sig- 
nificant uncertainties on these measurements, and for most 
sources any reasonable estimate of the error bars does not 
exclude a momentum flux that is close to L/c. Moreover, 
since the winds we observe now were launched some time 
ago, it is entirely possible that the present-day luminosity 
we measure for these sources is smaller than it was when the 
winds were launched. Given the observational uncertainties, 
we cannot conclude that giant clump winds cannot be radia- 
tively driven, only that, if they are, either their momentum 
fluxes must have been overestimated or their luminosities at 
the point of wind launching underestimated. 

It is interesting to ask how these conclusions compare 
with those of prior authors, and why they differ. Most no- 
tably, our conclusion that dust-reprocessed radiation is un- 
likely to be a significant factor in launching winds or disrupt- 
ing massive clusters is inconsistent with those of a number 



of authors, including Murray et al. (2010 2011), Hopkins 
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et al. (2011 2012), and Genel et al. (2012). We can under- 



stand the difference by examining how the radiation force 
imparted to the matter is computed in these models ver- 
sus in our simulations. Murray et al. ( 2010[ |2011 ) treat the 
matter as a one-dimensional thin shell, compute the result- 
ing gas temperature and thus the opacity, and compute the 
radiation force by multiplying the flux by that opacity. Simi- 
larly, |Hopkins_e^_aL ( 2011| 2012) adopt a uniform, high dust 



opacity (/^ = 5 cm g ) and compute the radiation force by 
multiplying that opacity by the total radiation flux. |Genel| 



et al. (2012) use a subgrid model that does not explicitly ac- 



count for radiation forces or any other feedback effect, but 
manually injects winds at a strength that is tuned to match 
the results of the [Murray et al. and [Hopkins et al.| models. 

Our numerical results show that these approaches can 
produce a large overestimate of the radiation force and the 
trapping factor, because they omit two critical effects re- 
vealed by our simulations. First, these models ignore the 
anti-correlation between radiation flux and gas density pro- 
duced by RRTI, which causes the flux seen by the bulk 
of the matter to be significantly smaller than the volume- 
averaged flux. Second, they miss the effect that, because 
the flux is low within the bulk of the matter, the gas and 
radiation temperatures are also low. This reduces the opac- 
ity within the bulk of the matter, further weakening the 
matter radiation-coupling. These two effects together yield 
values of /trap significantly smaller than suggested by previ- 
ous one-dimensional estimates. Given these results, it seems 



necessary to recompute the models of Murray et al.[ ([2010| 



2011), Hopkins et al. (2011 2012), and Genel et al.[ ([2012| 



using the approximate fitting formula for /trap that we have 
derived. 



5.4 Relation to Dusty Star Winds 

Although our work is focused on the problem of star cluster 
and galactic winds, the general problem of radiative driv- 
ing of dusty gas also arises in the context of winds from 
dusty late-type stars (e.g.' Goldreich Scoville|1976 Habing 
|199 6). Before proceeding, it is important to point out a sig- 
nificant way in which this problem differs from our work here 
and in Paper I. At the low temperatures typical of interstel- 
lar gas even in intens ely star-formin g galaxies, the opacity 
roughly oc (e.g. Sem enov et al.[200 3), as we use in our 
models. This means that the opacity general drops mono- 
tonically with height in an atmosphere. In the case of dusty 
stellar winds, on the other hand, where temperatures are 
near the grain sublimation temperature, the opacity is much 
more complex and non-monotonic, both due to grain forma- 
tion, destruction, and drift relative to the gas, and because 
even for a constant grain population the opacity varies are 



roughly (i.e. constant) rather than at temperatures 
close to the grain sublimation temperature. 

With this caveat aside, we note that one- dimensional 
models have been reasonably successful at reproducing many 



observations of massive star winds (e.g. Ivezic & Elitzur 



1995| |Ivezic k E litzur||2010| [Elitzur &: Ivezic[ [2001), impT^ 
ing that RRT instability may not be critical for these stars. 
However, the observationally-inf erred momenta of dusty star 
winds are usually below L/c (e.g. Groenewegen et al.|2009 ), 
and that even the highest inferred momenta are no more 
than ^ lOL/c (e.g. Elitzur & Ivezic 2001 ), which in turn sug- 



gests that these stars are not in the regime where we require 
strong amplification of the force by radiative trapping that 
might be inhibited by RRTI. We tentatively conclude that 
there is no contradiction between the observational results 
for AGB stars and our simulation results, but the problem is 
clearly worthy of further investigation, especially given the 
caveat above. 



6 SUMMARY 

In this paper we analyze the properties of optically thick 
radiation pressure- driven dusty winds. We consider the ide- 
alized problem of a column of material through which a 
specified radiation flux is passed. We first show that such 
a system is characterized by a single dimensionless num- 
ber, T*, the optical depth of the matter computed using 
the opacity at the dust photosphere, and that this param- 
eter will determine the rate at which the matter column 
absorbs momentum from the radiation field. We then use 
radiation-hydrodynamic simulations to measure this mo- 
mentum transfer rate. We find that, after one to a few dy- 
namical times, radiation Ray leigh- Taylor instability (RRTI) 
drives the gas into a configuration where most of the matter 
is in dense filaments aligned along the direction of the radi- 
ation flux, while most of the radiation flux passes through 
channels of reduced optical depth between the filaments. 
This configuration minimizes matter-radiation interaction, 
and thus limits the rate at which matter can take up mo- 
mentum from the radiation field. 

We combine this result with the result from Paper I, 
where we considered irradiated columns of matter confined 
by gravity, to produce a fitting formula for the behavior of 
irradiated, gravity- confined dusty gas layers. The behavior 
of these structures is determined by r* and by /e,*, the 
ratio of radiative and gravitational forces at the dust pho- 
tosphere. We identify three regimes of behavior depending 
on the values of these parameters. At a given r*, there is a 
critical value of /e,* below which radiation is dynamically 
unimportant. For values of /e,* above the critical value but 
below unity, RRTI sets in and makes the gas turbulent, but 
does not produce a wind. Only for /e,* > 1 is there a wind, 
and even in the limit where gravity provides negligible con- 
finement of that wind, the wind momentum flux is roughly 
1 + 0.5t* times the radiation momentum flux. 

We then consider the implications of these results for 
star-forming clusters and galaxies. For observed clusters and 
galaxies, our results suggest that dust-reprocessed radiation 
is unlikely to be able to drive winds and eject matter. The 
direct radiation field may still be able to launch winds, but 
only in systems where its momentum alone is sufficient to 
overcome gravity, without significant amplification by radia- 
tive trapping. 
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APPENDIX A: INITIAL CONDITIONS FOR 
RUN T30 

For run T30, we do not have initial conditions from Paper I 
because we did not perform any runs with r* =30. To gen- 
erate such conditions, we run a simulation with r* = 30 and 
/e,* = 0.05 (i.e. with gravity turned on) following the same 
procedure as for all other runs described in Paper I. We refer 
to this run as T30F0.05. As in the other runs from Paper 
I, we perform the simulation at a resolution Ax = 0.5/i*, in 
a computational domain of 1024 x 16384 cells, correspond- 
ing to a size of 512/i* x 8192/i*. We run the simulation for 
a time t = 75t*, by which point a turbulent flow is fully 
developed. To produce initial conditions for run T30 in this 
paper, we must rescale the results of run T30F0.05, because 
ha/h^ = 600, so that the resolution of run T30F0.05 is 
Ax = 8.3 X 10"^/ia. This is so high that it would be im- 
possible to advance the run for a time comparable to ta- 
We therefore downsample the output at the final time in 
run T30F0.05 by a factor of 32, producing a resolution of 
Ax = 0.027/ia- We also replicate the density, velocity, gas 
temperature, and radiation energy density fields 32 times 
in the horizontal direction; since run T30F0.05 has periodic 
boundary conditions, this is fully self-consistent. The result 
is a cube of initial conditions that is 1024 x 512 cells in size, 
at a resolution Ax = 0.027/ia, corresponding to a physical 
size 27.3/ia x 13.6/ia- We use this state as the initial con- 
dition for run T30, extending the computational domain in 
the vertical direction exactly as for the other runs described 
in Section [3l 
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